*-------------------------------------------------------------------------------
* bilevel b&b algorithm
*-------------------------------------------------------------------------------

set ii "leader problem variables' number"   /ii1/;
set i  "follower problem variables' number" /i1/;

parameter xp(ii);

variables
    f1, f2, f1r, f2r, f2p;

integer variables
    x(ii), xr(ii), y(i), yp(i);

positive variables
    yr(i);

x.lo(ii) = 0;
y.lo(i) = 0;
xr.lo(ii) = 0;
yp.lo(i) = 0;

equations
    obj1, obj2, e1, e2, e3,
    obj1r, obj2r, e1r, e2r, e3r,
    obj2f, e1f, e2f, e3f;

obj1..  f1 =e= -x('ii1') - 2*y('i1');
obj2..  f2 =e= y('i1');
e1..    -x('ii1') + 2.5*y('i1') =l= 3.75;
e2..    x('ii1') + 2.5*y('i1') =g= 3.75;
e3..    2.5*x('ii1') + y('i1') =l= 8.75;

obj1r.. f1r =e= -xr('ii1') - 2*yr('i1');
obj2r.. f2r =e= yr('i1');
e1r..   -xr('ii1') + 2.5*yr('i1') =l= 3.75;
e2r..   xr('ii1') + 2.5*yr('i1') =g= 3.75;
e3r..   2.5*xr('ii1') + yr('i1') =l= 8.75;

obj2f.. f2p =e= yp('i1');
e1f..   -xp('ii1') + 2.5*yp('i1') =l= 3.75;
e2f..   xp('ii1') + 2.5*yp('i1') =g= 3.75;
e3f..   2.5*xp('ii1') + yp('i1') =l= 8.75;

model high_point /obj1, e1, e2, e3/;

$echo bilevel x max f2r y obj2r, e1r, e2r, e3r > "%emp.info%"
model relaxed /obj1r, obj2r, e1r, e2r, e3r/;
option emp = lindo;

model follower /obj2f, e1f, e2f, e3f/;

*-------------------------------------------------------------------------------
* datastructure
*-------------------------------------------------------------------------------

set inode              "maximum size of integer node pool" /inode1*inode1000/;
parameter bound(inode) "inode n will have an obj <= bound(n)";
parameter iubound(inode,i)   "upper bound in this inode"; 
parameter ilbound(inode,i)   "lower bound in this inode";
parameter iiubound(inode,ii) "x upper bound in this inode"; 
parameter iilbound(inode,ii) "x lower bound in this inode";
set inewnode(inode)    "new inode";
set iwaiting(inode)    "waiting inode list";
set icurrent(inode)    "current inode";
parameter log(inode,*) "logging information";
scalar bestfound       "solution of follower problem (feasible)" /-INF/;
scalar bestpossible    "solution of high point problem (upper)"  /+INF/;
scalar done            "branch terminate"                        /0/;
scalar idone           "integer branch terminate"                /0/;
scalar objh            "objective of high_point problem";
scalar objc            "objective of relaxed problem";
scalar objf            "objective of follower's problem";

scalar first;
scalar first2;
scalar int;

parameter bn(i), bnn(ii);
scalar bnd, bndd;

scalar M /100/;

scalar flag /1/;

set w(inode);

alias(n,inode);

parameter inodenumber(inode);
inodenumber(inode) = ord(inode);

* add root inode to the waiting inode list

iwaiting('inode1') = yes;
icurrent('inode1') = yes;
inewnode('inode1') = yes;
bound('inode1') = +INF;

ilbound('inode1',i) = -M;
iubound('inode1',i) = M;
iilbound('inode1',ii) = -M;
iiubound('inode1',ii) = M;

loop(inode$(not idone),

* node selection
    bestpossible = smax(iwaiting(n), bound(n));
    icurrent(n) = no;
    icurrent(iwaiting(n))$(bound(n) = bestpossible) = yes;
    first = 1;

* only consider first in set current
    loop(icurrent$first,
        first = 0;


        log(inode,'inode') = inodenumber(icurrent);
        log(inode,'ub') = bestpossible;

* remove current inode from waiting list
        iwaiting(icurrent) = no;

* clear bounds
        x.lo(ii) = -M;
        x.up(ii) = M;
        xr.lo(ii) = -M;
        xr.lo(ii) = M;
        y.lo(i) = -M;
        y.up(i) = M;
        yr.lo(i) = -M;
        yr.up(i) = M;

* set appropriate bounds

        bn(i) = ilbound(icurrent,i) - iubound(icurrent,i);
        bnn(ii) = iilbound(icurrent,ii) - iiubound(icurrent,ii);
        bnd = smax(i,bn(i));
        bndd = smax(ii,bnn(ii));
        
        abort$(bnd gt 0 or bndd gt 0) "lower bound > upper bound";

        x.lo(ii) = iilbound(icurrent,ii);
        x.up(ii) = iiubound(icurrent,ii);
        xr.lo(ii) = iilbound(icurrent,ii);
        xr.up(ii) = iiubound(icurrent,ii);
        y.lo(i) = ilbound(icurrent,i);
        y.up(i) = iubound(icurrent,i);
        yr.lo(i) = ilbound(icurrent,i);
        yr.up(i) = iubound(icurrent,i);
        
        log(inode,'ylo') = yr.lo('i1');
        log(inode,'yup') = yr.up('i1');
        log(inode,'xup') = xr.up('ii1');
        log(inode,'xlo') = xr.lo('ii1');
        

* solve high point problem
        solve high_point using mip maximizing f1;

*        log(inode,'ssh') = high_point.solvestat;

* check for optimal solution
        abort$(high_point.solvestat <> 1) "Solver did not return ok for high point problem"
        if(high_point.modelstat eq 1 or high_point.modelstat eq 2,

* get objective
            objh = f1.l;
            bestpossible = objh;
            log(inode,'objh') = objh;

* solve relaxed problem
            solve relaxed using emp maximizing f1r;

*            log(inode,'ssr') = relaxed.solvestat;

            objc = f1r.l;
            log(inode,'objc') = objc;

* check for optimal solution
            abort$(relaxed.solvestat <> 1) "Solver did not return ok for relaxed problem";
            if(relaxed.modelstat eq 1 or relaxed.modelstat eq 2,

* check integrality
                int = smax(i, ceil(yr.l(i))-yr.l(i));
                if(int eq 0,
                    log(inode,'integer') = 1;

* fix leader's variables and solve follower's problem
                    xp(ii) = x.l(ii);
                    solve follower using mip maximizing f2p;
                    
                    abort$(follower.solvestat <> 1) "Solver did not return ok for follower problem";
                    if(follower.modelstat eq 1 or follower.modelstat eq 2,
* update bestfound
                        objf = f2p.l;

********************************************************************************
* change this bestfound for different problems
                        bestfound = -xp('ii1') - 2*yp.l('i1');
********************************************************************************

                        log(inode,'xp') = xp('ii1');
                        log(inode,'yp') = yp.l('i1');

                        log(inode,'bestfound') = bestfound;

* remove all waiting inodes with bestpossible < bestfound
                        w(n) = no; w(iwaiting) = yes;
                        iwaiting(w)$(bound(w) < bestfound) = no;

* integer branch
                        first2 = 1;

                        if(flag eq 1,
                            ilbound(icurrent,i) = -M;
                            iubound(icurrent,i) = M;
                        );

                        loop(ii$first2,
                            first2 = 0;
* create 2 new nodes for x
                            if(iiubound(icurrent,ii) gt xr.l(ii),
                                inewnode(n) = inewnode(n-1);
                                iubound(inewnode,i) = iubound(icurrent,i);
                                ilbound(inewnode,i) = ilbound(icurrent,i);
                                iiubound(inewnode,ii) = xr.l(ii) - 1;
                                iilbound(inewnode,ii) = iilbound(icurrent,ii);
                                bound(inewnode) = objh;
                                iwaiting(inewnode) = yes;
                            );

                            if(iilbound(icurrent,ii) lt xr.l(ii),
                                inewnode(n) = inewnode(n-1);
                                iubound(inewnode,i) = iubound(icurrent,i);
                                ilbound(inewnode,i) = ilbound(icurrent,i);
                                iiubound(inewnode,ii) = iiubound(icurrent,ii);
                                iilbound(inewnode,ii) = xr.l(ii) + 1;
                                bound(inewnode) = objh;
                                iwaiting(inewnode) = yes;
                            );
                        );

                        first2 = 1;

                        loop(i$first2,
                            first2 = 0;

* create 2 new nodes for y
                            if(iubound(icurrent,i) gt yr.l(i),
                                inewnode(n) = inewnode(n-1);
                                iubound(inewnode,i) = yr.l(i) - 1;
                                ilbound(inewnode,i) = ilbound(icurrent,i);
                                iiubound(inewnode,ii) = iiubound(icurrent,ii);
                                iilbound(inewnode,ii) = iilbound(icurrent,ii);
                                bound(inewnode) = objh;
                                iwaiting(inewnode) = yes;
                            );

                            if(ilbound(icurrent,i) lt yr.l(i),
                                inewnode(n) = inewnode(n-1);
                                iubound(inewnode,i) = iubound(icurrent,i);
                                ilbound(inewnode,i) = yr.l(i) + 1;
                                iiubound(inewnode,ii) = iiubound(icurrent,ii);
                                iilbound(inewnode,ii) = iilbound(icurrent,ii);
                                bound(inewnode) = objh;
                                iwaiting(inewnode) = yes;
                            );
                        );
                    );

                    flag = 0;

                else
                    log(inode,'integer') = -1;

* branch
                    first2 = 1;

                    if(flag eq 0,
                        ilbound(icurrent,i) = -M;
                        iubound(icurrent,i) = M;
                    );

                    loop(i$first2,
                        first2 = 0;
* create 2 new nodes
                        if(iubound(icurrent,i) gt ilbound(icurrent,i),
                            inewnode(n) = inewnode(n-1);
                            iubound(inewnode,i) = ceil(yr.l(i)) - 1;
                            ilbound(inewnode,i) = ilbound(icurrent,i);
                            iiubound(inewnode,ii) = iiubound(icurrent,ii);
                            iilbound(inewnode,ii) = iilbound(icurrent,ii);
                            bound(inewnode) = objh;
                            iwaiting(inewnode) = yes;
                        );
                        
                        if(iubound(icurrent,i) gt ilbound(icurrent,i),
                            inewnode(n) = inewnode(n-1);
                            iubound(inewnode,i) = iubound(icurrent,i);
                            ilbound(inewnode,i) = ceil(yr.l(i));
                            iiubound(inewnode,ii) = iiubound(icurrent,ii);
                            iilbound(inewnode,ii) = iilbound(icurrent,ii);
                            bound(inewnode) = objh;
                            iwaiting(inewnode) = yes;
                        );
                    );

                    flag = 1;
                );
            );
        );
    );

    
    
* are there new waiting inode?
    idone$(card(iwaiting) = 0) = 1;

* display log
    display log;
);
